Iris Pretzlaff at Hamburg University studied hibernation of dormice with the goal of understanding how global warming will affect physiology during hibernation and how this can affect population viability
Hypotheses:
Iris observed 12 and 13 individuals over winters of 2009 and 2011 in wooden boxes, kept outside in their natural conditions. Skin temperature was measured with data loggers. Torpor and arousal bouts were assigned based on skin temperature: Torpor was defined as any time when skin temperature decreased < 30°C (Tskin of active, normothermic animals is around 36°C) and an arousal when skin temperature increased > 30°C between torpor bouts, or when the difference Tskin-Tnest was > 30°C.
Metabolic rate was measured as oxygen consumption (ml \(O_{2}\)/h) using open-flow respirometry. For more details see Pretzlaff et al. (2021).
We focus on testing the hypothesis about the probability of arousal.
xyplot(state ~ Tamb|ID, data = df)
xyplot(state ~ HoursBegin|ID, data = df)
s2009 <- df[df$Year %in% c(2009, 2010), ]
s2011 <- df[df$Year %in% c(2011, 2012), ]
xyplot(state ~ idate | ID, data = s2009)
xyplot(state ~ idate | ID, data = s2011)
We rescale the data for better model convergence. Feel free to try and fit the model without this step (I already did that) and see what happens.
df <- df[! (is.na(df$HoursBegin) | is.na(df$Tamb)), ]
min(df$HoursBegin)
## [1] 19
# coding year as a factor
df$Y <- 2009
df$Y[df$Year %in% c(2011, 2012)] <- 2011
df$Y <- as.factor(df$Y)
### rescaling time since hibernation and temperature for better convergence
df$HoursScale <- (df$HoursBegin - mean(df$HoursBegin, na.rm = T)) / sd(df$HoursBegin, na.rm = T)
df$HoursScale2 <- df$HoursScale^2
hist(df$HoursScale,
xlab = 'Standarized time since begin hibernation',
main = NULL, col = 'black')
df$T_scale <- (df$Tamb - mean(df$Tamb, na.rm = T)) / sd(df$Tamb, na.rm = T)
df$Tscale2 <- df$T_scale ^ 2
hist(df$T_scale, xlab = 'Standardized temperature',
main = NULL, col ='black')
mod0_probA <- glm(state ~ HoursScale + HoursScale2 +
T_scale + Tscale2 + Y,
family = binomial(link = 'logit'),
data = df)
summary(mod0_probA)
##
## Call:
## glm(formula = state ~ HoursScale + HoursScale2 + T_scale + Tscale2 +
## Y, family = binomial(link = "logit"), data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9353 -0.5560 -0.3860 -0.2337 2.7603
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.35664 0.05807 -57.803 < 2e-16 ***
## HoursScale 0.01895 0.02040 0.929 0.353
## HoursScale2 0.72028 0.02783 25.883 < 2e-16 ***
## T_scale 0.45224 0.03498 12.930 < 2e-16 ***
## Tscale2 0.11389 0.01913 5.954 2.61e-09 ***
## Y2011 0.95357 0.04988 19.119 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 14940 on 16579 degrees of freedom
## Residual deviance: 11951 on 16574 degrees of freedom
## AIC: 11963
##
## Number of Fisher Scoring iterations: 6
Any idea how we can conduct model
diagnostics for this model?
What about calculating dispersion, do you
remember how we can do a quick check whether overdispersion is
potentially an issue?
Here we assume that time since the beginning of hibernation and temperature have same effects across all individuals, but that there is some random variation among individuals in the intercept for i) time since begin of hibernation and ii) temperature for each individual. We use glmmTMB package to fit these models because we want to be able to plot the predictions in the end (+ ideally fit autocorrelation per individual)
modRandID_prob <- glmmTMB(state ~ HoursScale + HoursScale2 +
T_scale + Tscale2 +
Y + (1|ID),
family = binomial,
data=df)
summary(modRandID_prob)
## Family: binomial ( logit )
## Formula:
## state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + Y + (1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 11147.4 11201.4 -5566.7 11133.4 16573
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.6676 0.8171
## Number of obs: 16580, groups: ID, 25
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.55700 0.23680 -15.021 < 2e-16 ***
## HoursScale 0.11709 0.02278 5.141 2.74e-07 ***
## HoursScale2 0.80785 0.03111 25.969 < 2e-16 ***
## T_scale 0.48802 0.03645 13.387 < 2e-16 ***
## Tscale2 0.14762 0.01906 7.746 9.49e-15 ***
## Y2011 0.92287 0.33263 2.774 0.00553 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Diagnostics of the fitted random-intercept model. Let us use DHARMa package.
sim <- simulateResiduals(modRandID_prob, plot = T)
testUniformity(sim, plot = F)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.011804, p-value = 0.0197
## alternative hypothesis: two-sided
testQuantiles(sim, plot = T)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value < 2.2e-16
## alternative hypothesis: both
testDispersion(sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.92247, p-value = 0.296
## alternative hypothesis: two.sided
# plot residuals vs our predictors
plotResiduals(sim, form = df$T_scale)
plotResiduals(sim, form = df$HoursScale)
What does this diagnostics plots tell us?
Test for temporal autocorrelation since the measurements were taken over time.
sim_ID10 <- recalculateResiduals(sim, sel = df$ID == 10)
testTemporalAutocorrelation(sim_ID10, time = df$Timecon[df$ID == '10'])
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.6427, p-value = 1.239e-07
## alternative hypothesis: true autocorrelation is not 0
Let us use a similar function we used in the lecture to test temporal autocorrelation for each individual in the dataset.
levels(df$ID)
## [1] "10" "11A" "11B" "12" "14A" "14B" "2" "21" "22" "2222"
## [11] "23" "24A" "24B" "38A" "38B" "38C" "4" "57" "64" "70"
## [21] "79" "86" "97A" "97B" "98"
fun_autocor(sim, df, ID = '10')
## [[1]]
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.6427, p-value = 1.239e-07
## alternative hypothesis: true autocorrelation is not 0
##
##
## [[2]]
## DurbW p ID
## DW 1.64269 1.23855e-07 10
dat_autoc <- data.frame(DurbW =numeric(), p = numeric())
dat_autoc <- rbind(dat_autoc, fun_autocor(sim, df, 10)[[2]])
for(i in unique(df$ID)[-1]){
dat_autoc <- rbind(dat_autoc, fun_autocor(sim, df, ID = i)[[2]])
}
dat_autoc$sign <- ifelse(dat_autoc$p < 0.05,'Yes','No')
sum(dat_autoc$sign == 'Yes')
## [1] 25
length(levels(df$ID))
## [1] 25
So, autocorrelation is an issue
Attention: we assume we have time running continuously for each individual. But in fact, for some of them there are little gaps in the records (for temperature and metabolic rate)
df$Time_fac <- as.factor(df$Timecon)
modRandAR_prob <- glmmTMB(state ~ HoursScale + HoursScale2 +
T_scale + Tscale2 +
Y + (1|ID) + ar1(Time_fac - 1|ID),
family = binomial,
data=df)
summary(modRandAR_prob)
## Family: binomial ( logit )
## Formula:
## state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + Y + (1 |
## ID) + ar1(Time_fac - 1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 3893.1 3962.6 -1937.6 3875.1 16571
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 8.467e-07 9.201e-04
## ID.1 Time_fac1 1.937e+03 4.401e+01 0.93 (ar1)
## Number of obs: 16580, groups: ID, 25
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.3817 6.3070 -9.574 < 2e-16 ***
## HoursScale 1.1673 0.9665 1.208 0.22710
## HoursScale2 13.7559 1.6632 8.271 < 2e-16 ***
## T_scale 8.3990 1.3541 6.203 5.55e-10 ***
## Tscale2 1.6549 0.5481 3.020 0.00253 **
## Y2011 19.2246 3.1915 6.024 1.71e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
How do you interpret these results? What do you think about this model compared to the one without autocorrelation?
simAR1 <- simulateResiduals(modRandAR_prob, plot = T)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
testUniformity(simAR1, plot = F)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.066222, p-value < 2.2e-16
## alternative hypothesis: two-sided
testQuantiles(simAR1, plot = T)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value < 2.2e-16
## alternative hypothesis: both
testDispersion(simAR1)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.6603, p-value < 2.2e-16
## alternative hypothesis: two.sided
# plot residuals vs our predictors
plotResiduals(simAR1, form = df$T_scale)
plotResiduals(simAR1, form = df$HoursScale)
What are main issues with model diagnostics? What can be done to fit a model that would match the data structure better?
We have fitted the models assuming the slopes are exactly the same across individuals and only intercepts differ. How realistic is that? Should we have tried a different random structure?
From the note on DHARMa documentation:
Standard DHARMa simulations from models with (temporal / spatial / phylogenetic) conditional autoregressive terms will still have the respective temporal / spatial / phylogenetic correlation in the DHARMa residuals, unless the package you are using is modelling the autoregressive terms as explicit REs and is able to simulate conditional on the fitted REs. This has two consequences: 1. If you check the residuals for such a model, they will still show significant autocorrelation, even if the model fully accounts for this structure. 2. Because the DHARMa residuals for such a model are not statistically independent any more, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.
Attention! This is very slow and takes long time to run. I suggest you run it after the session
##DO not run
simAR1_rot <- simulateResiduals(modRandAR_prob, plot = T, rotation = 'estimated')
testUniformity(simAR1_rot, plot = F)
testQuantiles(simAR1_rot, plot = T)
testDispersion(simAR1_rot)
# plot residuals vs our predictors
plotResiduals(simAR1_rot, form = df$T_scale)
plotResiduals(simAR1_rot, form = df$HoursScale)
Conclusion: know what function is doing exactly. Do not just copy-paste the codes!
To test the hypothesis that both temperature and time have a quadratic effect on probability of arousal, we fit the models without each quadratic term in turn, and run LRTs.
modRandAR_no2T_prob <- glmmTMB(state ~ HoursScale + HoursScale2 +
T_scale +
Y + (1 | ID) + ar1(Time_fac - 1|ID),
family = binomial,
data=df)
## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence
## (7). See vignette('troubleshooting')
modRandAR_no2Hour_prob <- glmmTMB(state ~ HoursScale +
T_scale + Tscale2 +
Y + (1 | ID) + ar1(Time_fac - 1|ID),
family = binomial,
data=df)
## Warning in fitTMB(TMBStruc): Model convergence problem; false convergence (8).
## See vignette('troubleshooting')
anova(modRandAR_prob, modRandAR_no2T_prob) ## significant quadratic effect
## Data: df
## Models:
## modRandAR_no2T_prob: state ~ HoursScale + HoursScale2 + T_scale + Y + (1 | ID) + ar1(Time_fac - , zi=~0, disp=~1
## modRandAR_no2T_prob: 1 | ID), zi=~0, disp=~1
## modRandAR_prob: state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + Y + (1 | , zi=~0, disp=~1
## modRandAR_prob: ID) + ar1(Time_fac - 1 | ID), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modRandAR_no2T_prob 8 3905.8 3967.5 -1944.9 3889.8
## modRandAR_prob 9 3893.1 3962.6 -1937.6 3875.1 14.65 1 0.0001294
##
## modRandAR_no2T_prob
## modRandAR_prob ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note we have ‘Warning about false convergence’ for the model without the quadratic effect of time. Inspect its output and think what are the possible causes and what can be done?
summary(modRandAR_no2Hour_prob)
## Family: binomial ( logit )
## Formula:
## state ~ HoursScale + T_scale + Tscale2 + Y + (1 | ID) + ar1(Time_fac -
## 1 | ID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 3997.9 4059.7 -1991.0 3981.9 16572
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.661e-06 0.001289
## ID.1 Time_fac1 1.242e+03 35.242402 0.93 (ar1)
## Number of obs: 16580, groups: ID, 25
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -35.4141 3.8188 -9.274 < 2e-16 ***
## HoursScale -2.9674 0.9732 -3.049 0.00229 **
## T_scale 11.2666 1.2090 9.319 < 2e-16 ***
## Tscale2 2.2701 0.3733 6.081 1.19e-09 ***
## Y2011 12.4581 2.7251 4.572 4.84e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Prepare a new data frame. First predict on the scale of the linear predictor.
min(df$T_scale)
## [1] -4.340506
min(df$HoursScale)
## [1] -2.060011
new_dat <- expand.grid('ID' = unique(df$ID),
'Y' = unique(df$Y),
T_scale = seq(from = min(df$T_scale), to =
max(df$T_scale), length.out = 20),
HoursScale = seq(from = min(df$HoursScale), to =
max(df$HoursScale), length.out = 20))
new_dat$HoursScale2 <- new_dat$HoursScale^2
new_dat$Tscale2 <- new_dat$T_scale^2
head(new_dat)
## ID Y T_scale HoursScale HoursScale2 Tscale2
## 1 10 2011 -4.340506 -2.060011 4.243644 18.83999
## 2 11A 2011 -4.340506 -2.060011 4.243644 18.83999
## 3 11B 2011 -4.340506 -2.060011 4.243644 18.83999
## 4 12 2011 -4.340506 -2.060011 4.243644 18.83999
## 5 14A 2011 -4.340506 -2.060011 4.243644 18.83999
## 6 14B 2011 -4.340506 -2.060011 4.243644 18.83999
new_dat$Time_fac <- as.factor(rep('1', nrow(new_dat))) ## keeping just to one time point, to avoid variaiton here
new_dat$pred_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'link')
new_dat <- new_dat[order(new_dat$ID, new_dat$Y, new_dat$T_scale), ]
ggplot(data = new_dat, aes(x = T_scale, y = pred_fixTMB, group = as.factor(HoursScale), col = as.factor(HoursScale))) +
geom_point() + theme_bw() +
facet_grid(cols = vars(Y)) +
xlab('Standardized temperature') + ylab('Linear predictor')
ggplot(data = new_dat, aes(x = HoursScale, y = pred_fixTMB, group = as.factor(HoursScale), col = as.factor(T_scale))) +
geom_point() + theme_bw() +
facet_grid(cols = vars(Y)) +
xlab('Standardized time sing begin of hibernation') + ylab('Linear predictor')
What is the probability to arouse close to the end of hibernation
max(df$HoursScale)
## [1] 1.841177
new_dat <- expand.grid('ID' = unique(df$ID),
'Y' = unique(df$Y),
T_scale = seq(from = min(df$T_scale), to =
max(df$T_scale), length.out = 20),
HoursScale = 1.6) ## set time to median, to avoid additional variaiton
new_dat$HoursScale2 <- new_dat$HoursScale^2
new_dat$Tscale2 <- new_dat$T_scale^2
head(new_dat)
## ID Y T_scale HoursScale HoursScale2 Tscale2
## 1 10 2011 -4.340506 1.6 2.56 18.83999
## 2 11A 2011 -4.340506 1.6 2.56 18.83999
## 3 11B 2011 -4.340506 1.6 2.56 18.83999
## 4 12 2011 -4.340506 1.6 2.56 18.83999
## 5 14A 2011 -4.340506 1.6 2.56 18.83999
## 6 14B 2011 -4.340506 1.6 2.56 18.83999
new_dat$Time_fac <- as.factor(rep('1', nrow(new_dat))) ## keeping just to one time point, to avoid variaiton here
new_dat$pred_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'response')
ggplot(data = new_dat, aes(x = T_scale, y = pred_fixTMB)) +
geom_line() + theme_bw() +
facet_grid(cols = vars(Y)) +
xlab('Standardized temperature') + ylab('Probability of arousal')
Add CI to the prediction on probability to arouse from hibernation
new_dat$pred_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'response')
new_dat$pre_SE_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, se.fit = TRUE, type = 'response')$se.fit
new_dat$pred_upper_fix <- new_dat$pred_fixTMB + qnorm(0.975)*new_dat$pre_SE_fixTMB
new_dat$pred_lower_fix <- new_dat$pred_fixTMB - qnorm(0.975)*new_dat$pre_SE_fixTMB
head(new_dat)
## ID Y T_scale HoursScale HoursScale2 Tscale2 Time_fac pred_fixTMB
## 1 10 2011 -4.340506 1.6 2.56 18.83999 1 8.682929e-05
## 2 11A 2011 -4.340506 1.6 2.56 18.83999 1 8.682929e-05
## 3 11B 2011 -4.340506 1.6 2.56 18.83999 1 8.682929e-05
## 4 12 2011 -4.340506 1.6 2.56 18.83999 1 8.682929e-05
## 5 14A 2011 -4.340506 1.6 2.56 18.83999 1 8.682929e-05
## 6 14B 2011 -4.340506 1.6 2.56 18.83999 1 8.682929e-05
## pre_SE_fixTMB pred_upper_fix pred_lower_fix
## 1 0.0007073774 0.001473263 -0.001299605
## 2 0.0007073774 0.001473263 -0.001299605
## 3 0.0007073774 0.001473263 -0.001299605
## 4 0.0007073774 0.001473263 -0.001299605
## 5 0.0007073774 0.001473263 -0.001299605
## 6 0.0007073774 0.001473263 -0.001299605
# predictions for the mean response across species
ggplot(data = new_dat, aes(x = T_scale, y = pred_fixTMB)) +
theme_bw() + facet_grid(cols = vars(Y)) +
geom_ribbon(data = new_dat, aes(x = T_scale, ymin = pred_lower_fix, ymax = pred_upper_fix), alpha = 0.1) +
geom_line(lty = 3) +
xlab('Standardized temperature') + ylab('Probability of arousal')
Attention: we only used predictions with the fixed part of the model,
not including random effects at all!
Repetition: how do you go from the scale of the response to the linear predictor?
new_dat$pred_fixTMB_link <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'link')
new_dat$pred_trans <- exp(new_dat$pred_fixTMB_link)/(1 + exp(new_dat$pred_fixTMB_link))
ggplot(data = new_dat, aes(x = T_scale, y = pred_trans)) +
geom_line() + theme_bw() +
facet_grid(cols = vars(Y)) +
xlab('Standardized temperature') + ylab('Probability of arousal')